{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.14",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "kaggle": {
   "accelerator": "nvidiaTeslaT4",
   "dataSources": [],
   "dockerImageVersionId": 30787,
   "isInternetEnabled": true,
   "language": "python",
   "sourceType": "notebook",
   "isGpuEnabled": true
  }
 },
 "nbformat_minor": 5,
 "nbformat": 4,
 "cells": [
  {
   "id": "af4b048c0e49992d",
   "cell_type": "code",
   "source": "import torch\nfrom collections import OrderedDict\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.io\n#from pyDOE import lhs#拉丁超立方抽样\nfrom torch import nn\nimport time\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=UserWarning)  # 只忽略UserWarning类型的警告",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T07:12:44.360332Z",
     "iopub.execute_input": "2024-10-22T07:12:44.360725Z",
     "iopub.status.idle": "2024-10-22T07:12:44.366658Z",
     "shell.execute_reply.started": "2024-10-22T07:12:44.360686Z",
     "shell.execute_reply": "2024-10-22T07:12:44.365544Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:40.677214Z",
     "start_time": "2024-10-22T06:51:40.673203Z"
    }
   },
   "execution_count": 51,
   "outputs": []
  },
  {
   "id": "9e2f9d35aa33fbb9",
   "cell_type": "code",
   "source": "if torch.cuda.is_available():\n    device = torch.device('cuda')\nelse:\n    device = torch.device('cpu')\ndef setup_seed(seed):\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed_all(seed)\n    np.random.seed(seed)\n    #  random.seed(seed)\n    torch.backends.cudnn.deterministic = True\n\nsetup_seed(124)\ndtype = torch.double",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T07:12:44.582623Z",
     "iopub.execute_input": "2024-10-22T07:12:44.582997Z",
     "iopub.status.idle": "2024-10-22T07:12:44.595530Z",
     "shell.execute_reply.started": "2024-10-22T07:12:44.582959Z",
     "shell.execute_reply": "2024-10-22T07:12:44.594730Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:40.926021Z",
     "start_time": "2024-10-22T06:51:40.916997Z"
    }
   },
   "execution_count": 52,
   "outputs": []
  },
  {
   "id": "e501d8818e912d51",
   "cell_type": "code",
   "source": "#求导\ndef gradients(outputs, inputs):\n    return torch.autograd.grad(outputs, inputs, grad_outputs=torch.ones_like(outputs), create_graph=True)\n#类型转换\ndef to_numpy(input):\n    if isinstance(input, torch.Tensor):\n        return input.detach().cpu().numpy()\n    elif isinstance(input, np.ndarray):\n        return input\n",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T07:12:45.020677Z",
     "iopub.execute_input": "2024-10-22T07:12:45.021071Z",
     "iopub.status.idle": "2024-10-22T07:12:45.027126Z",
     "shell.execute_reply.started": "2024-10-22T07:12:45.021034Z",
     "shell.execute_reply": "2024-10-22T07:12:45.026190Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:41.142575Z",
     "start_time": "2024-10-22T06:51:41.126523Z"
    }
   },
   "execution_count": 53,
   "outputs": []
  },
  {
   "id": "79e6a971aae89fdf",
   "cell_type": "code",
   "source": "#构建网络layers\ndef Dnnlayers(input_layers,output_layers,hidden_layers,neural):\n    layers = []\n    for i in range(hidden_layers + 2):\n        if i == 0:\n            layers.append(input_layers)\n        elif i == hidden_layers + 1:\n            layers.append(output_layers)\n        else:\n            layers.append(neural)\n    return layers\n#初值\ndef IC(x, crhoL, cuL, cpL, crhoR, cuR, cpR):\n    N = x.shape[0]\n    rho_init = np.zeros((x.shape[0], 1))\n    u_init = np.zeros((x.shape[0], 1))\n    p_init = np.zeros((x.shape[0], 1))\n\n    # rho, p - initial condition\n    for i in range(N):\n        if (x[i, 1] <= 0.0):\n            rho_init[i] = crhoL\n            u_init[i] = cuL\n            p_init[i] = cpL\n        else:\n            rho_init[i] = crhoR\n            u_init[i] = cuR\n            p_init[i] = cpR\n    U_ic = np.hstack([rho_init, u_init, p_init])\n\n    return U_ic\n#外推边界\ndef BC(x, crhoL, cuL, cpL, crhoR, cuR, cpR):\n    N = x.shape[0]\n    rho_bc = np.zeros((x.shape[0], 1))\n    u_bc = np.zeros((x.shape[0], 1))\n    p_bc = np.zeros((x.shape[0], 1))\n\n    # rho, p - initial condition\n    for i in range(N):\n        if (x[i, 1] <= 0.0):\n            rho_bc[i] = crhoL\n            u_bc[i] = cuL\n            p_bc[i] = cpL\n        else:\n            rho_bc[i] = crhoR\n            u_bc[i] = cuR\n            p_bc[i] = cpR\n    U_bc = np.hstack([rho_bc, u_bc, p_bc])\n\n    return U_bc\n",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T07:12:45.280709Z",
     "iopub.execute_input": "2024-10-22T07:12:45.281086Z",
     "iopub.status.idle": "2024-10-22T07:12:45.293293Z",
     "shell.execute_reply.started": "2024-10-22T07:12:45.281048Z",
     "shell.execute_reply": "2024-10-22T07:12:45.292325Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:41.360999Z",
     "start_time": "2024-10-22T06:51:41.344299Z"
    }
   },
   "execution_count": 54,
   "outputs": []
  },
  {
   "id": "2b5f4f99b1b8181f",
   "cell_type": "code",
   "source": "class DNN(nn.Module):\n    def __init__(self, layers):\n        #继承父类\n        super(DNN, self).__init__()\n\n        #depth of network\n        self.depth = len(layers) - 1\n        #print(self.depth)\n\n        #activation of network\n        self.activation = nn.Tanh()\n\n        #create the neural network\n        layers_list = list()\n        for i in range(self.depth - 1):\n            layers_list.append(\n                ('layer_%d' % i, nn.Linear(layers[i], layers[i+1]))#create each network\n            )\n            layers_list.append(\n                ('activation_%d' % i, self.activation)  #create each activation of network\n            )\n        layers_list.append(\n            ('layer_%d' % (self.depth - 1), nn.Linear(layers[-2], layers[-1]))\n            #last network do not have activation\n        )\n        #创建一个有序字典，其中包含了从 layers_list 中获得的键值对，这在需要保持元素顺序的场景（如神经网络层的顺序）中非常有用。\n        layerDict = OrderedDict(layers_list)\n\n        #deploy layers\n        self.layers = nn.Sequential(layerDict)\n        #print(self.layers)\n\n    #forword network,output the result of network\n    def forword(self, x):\n        out = self.layers(x)\n        return out",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T07:12:45.561377Z",
     "iopub.execute_input": "2024-10-22T07:12:45.561759Z",
     "iopub.status.idle": "2024-10-22T07:12:45.570558Z",
     "shell.execute_reply.started": "2024-10-22T07:12:45.561720Z",
     "shell.execute_reply": "2024-10-22T07:12:45.569643Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:41.563292Z",
     "start_time": "2024-10-22T06:51:41.548006Z"
    }
   },
   "execution_count": 55,
   "outputs": []
  },
  {
   "id": "5a835b6157600a38",
   "cell_type": "code",
   "source": "class Pinns():\n    def __init__(self, layers):\n        self.dnn = DNN(layers).double().to(device)    \n    #forword network,output the result of network\n    def forword(self, x):\n        y = self.dnn.forword(x)\n        rho, u, p = y[:, 0:1], y[:, 1:2], y[:, 2:]\n        return rho, u, p\n    #function loss,残差点训练\n    def loss_pde(self, x_int):\n        y = self.dnn.forword(x_int)\n        rho, u, p = y[:, 0:1], y[:, 1:2], y[:, 2:]\n        #gamma = self.gamma\n\n        U2 = rho * u\n        U3 = 0.5 * rho * u ** 2 + p / (1.4 - 1.0)#0.4 #gamma - 1 = 0.4\n        # F1 = U2\n        F2 = rho * u ** 2 + p\n        F3 = u * (U3 + p)\n\n        # Gradients and partial derivatives\n        drho_g = gradients(rho, x_int)[0]\n        rho_t, rho_x = drho_g[:, :1], drho_g[:, 1:]\n\n        dU2_g = gradients(U2, x_int)[0]\n        U2_t, U2_x = dU2_g[:, :1], dU2_g[:, 1:]\n        dU3_g = gradients(U3, x_int)[0]\n        U3_t, U3_x = dU3_g[:, :1], dU3_g[:, 1:]\n        dF2_g = gradients(F2, x_int)[0]\n        F2_t, F2_x = dF2_g[:, :1], dF2_g[:, 1:]\n        dF3_g = gradients(F3, x_int)[0]\n        F3_t, F3_x = dF3_g[:, :1], dF3_g[:, 1:]\n        \n        f = (((rho_t + U2_x) ) ** 2).mean() + \\\n            (((U2_t + F2_x) ) ** 2).mean() + \\\n            (((U3_t + F3_x) ) ** 2).mean()  #\n        \n        return f\n    #initial loss function\n    def loss_ic(self, x_ic, U_ic):\n        y_ic = self.dnn.forword(x_ic)\n        rho_ic_nn, u_ic_nn, p_ic_nn = y_ic[:, 0:1], y_ic[:, 1:2], y_ic[:, 2:]\n        loss_ics = ((rho_ic_nn - U_ic[:,0:1]) ** 2).mean()+\\\n                   ((u_ic_nn - U_ic[:,1:2]) ** 2).mean()+\\\n                   ((p_ic_nn - U_ic[:,2:]) ** 2).mean()\n        return loss_ics\n    #boundary condition\n    def loss_bc(self, x_bc, U_bc):\n        y_bc = self.dnn.forword(x_bc)\n        rho_bc_nn, u_bc_nn, p_bc_nn = y_bc[:, 0:1], y_bc[:, 1:2], y_bc[:, 2:]\n        loss_bcs = ((rho_bc_nn - U_bc[:, 0:1]) ** 2).mean() + \\\n                   ((u_bc_nn - U_bc[:, 1:2]) ** 2).mean() + \\\n                   ((p_bc_nn - U_bc[:, 2:]) ** 2).mean()\n        return loss_bcs",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T07:12:45.888715Z",
     "iopub.execute_input": "2024-10-22T07:12:45.889121Z",
     "iopub.status.idle": "2024-10-22T07:12:45.905453Z",
     "shell.execute_reply.started": "2024-10-22T07:12:45.889082Z",
     "shell.execute_reply": "2024-10-22T07:12:45.904455Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:41.749721Z",
     "start_time": "2024-10-22T06:51:41.722083Z"
    }
   },
   "execution_count": 56,
   "outputs": []
  },
  {
   "id": "6be409ae6f6ebf92",
   "cell_type": "code",
   "source": "class conbine():\n\n    def __init__(self, layers1, layers2, layers3):\n\n        self.model1 = Pinns(layers1)\n\n        self.model2 = Pinns(layers2)\n\n        self.model3 = DNN(layers3).double().to(device)    \n\n    def loss_pde(self, x_int):   \n\n        index1 = torch.where(x_int[:, 1]<=x_predict[:,0])[0]\n\n        index2 = torch.where(x_int[:, 1]>x_predict[:,0])[0]\n\n        x1_int = x_int[index1]\n\n        x2_int = x_int[index2]\n\n        #loss_pde(self, x_int):\n\n        loss_pde1 = self.model1.loss_pde(x1_int)\n\n        loss_pde2 = self.model2.loss_pde(x2_int)\n\n        loss_pde = loss_pde1 + loss_pde2\n\n        return loss_pde\n\n    #loss_ic(self, x_ic, U_ic):\n\n    def loss_ic(self, x_ic, U_ic):        \n\n        index1 = torch.where(x_ic[:, 1]<=x_predict_ic[:,0])[0]\n\n        index2 = torch.where(x_ic[:, 1]>x_predict_ic[:,0])[0]\n\n        x1_ic = x_ic[index1]\n\n        x2_ic = x_ic[index2]\n\n        U1_ic = U_ic[index1]\n\n        U2_ic = U_ic[index2]\n\n        loss_ic1 = self.model1.loss_ic(x1_ic, U1_ic)\n\n        loss_ic2 = self.model2.loss_ic(x2_ic, U2_ic)\n\n        loss_ic = loss_ic1 + loss_ic2\n\n        return loss_ic\n\n    def loss_xic(self, x_ic):\n\n        loss_ic = ((x_predict_ic - x_ic[:, 0:1]) ** 2).mean()\n\n        return loss_ic\n\n    #loss_bc(self, x_bc, U_bc):\n\n    #R-H，D=u关系\n\n    def loss_rh(self, t_rh_int, x_start, x_end):\n\n        x_predict = model.model3.forword(t_rh_int)[:, 0:1]\n\n        index = torch.where((x_predict[:, 0]<=x_end)&(x_predict[:, 0]>=x_start))\n\n        dx = gradients(x_predict, t_rh_int)[0]\n\n        dx_t, dx_x = dx[:, 0:1], dx[:, 1:2]\n\n        dx_t = dx_t[index]\n\n        #\n\n        t_rh_int = t_rh_int[index]\n\n        rho1, u1, p1 = self.model1.forword(t_rh_int)\n\n        rho2, u2, p2 = self.model2.forword(t_rh_int)\n\n        loss_rh = ((u1 - u2) ** 2).mean()+ ((p1 - p2) ** 2).mean()+((dx_t - u1) ** 2).mean()+((dx_t - u2) ** 2).mean()+((dx_x) ** 2).mean()\n\n        return loss_rh\n\n    def loss_x(self):\n\n        x_int = t_rh_int_train[:,1:2]\n\n        loss_xy = ((x_predict-x_int)**2).mean()\n\n        return loss_xy\n\n       ",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:43.776167Z",
     "iopub.execute_input": "2024-10-22T08:06:43.776939Z",
     "iopub.status.idle": "2024-10-22T08:06:43.794732Z",
     "shell.execute_reply.started": "2024-10-22T08:06:43.776890Z",
     "shell.execute_reply": "2024-10-22T08:06:43.793846Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:41.921829Z",
     "start_time": "2024-10-22T06:51:41.904253Z"
    }
   },
   "execution_count": 119,
   "outputs": []
  },
  {
   "id": "84aca8ea777720d8",
   "cell_type": "code",
   "source": "#网格构建\nnum_x = 401\nnum_t = 201\nTstart = 0.0\nTend = 0.5\nXstart = -1.0\nXend = 1.0\nx = np.linspace(Xstart, Xend, num_x)\nt = np.linspace(Tstart, Tend, num_t)\nt_grid, x_grid = np.meshgrid(t, x)\nT = t_grid.flatten()[:, None]\nX = x_grid.flatten()[:, None]\nX_star = np.hstack((T.flatten()[:, None], X.flatten()[:, None]))\n# 时空边界\nlb = X_star.min(0)  # 时空下界，\nub = X_star.max(0)  # 时空上界\nprint('lb is :',lb, ', ub is :', ub)",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:44.310893Z",
     "iopub.execute_input": "2024-10-22T08:06:44.311561Z",
     "iopub.status.idle": "2024-10-22T08:06:44.326441Z",
     "shell.execute_reply.started": "2024-10-22T08:06:44.311521Z",
     "shell.execute_reply": "2024-10-22T08:06:44.325468Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:42.108302Z",
     "start_time": "2024-10-22T06:51:42.083245Z"
    }
   },
   "execution_count": 120,
   "outputs": []
  },
  {
   "id": "bb1f3322bfae61ab",
   "cell_type": "code",
   "source": "#初值\n\n#loss_ic(self, x_ic, U_ic):\n\nnum_i_train = 201\n\nid_ic = np.random.choice(num_x, num_i_train, replace=False)\n\nirhoL, iuL, ipL = 1.0, 1.0, 0.125\n\nirhoR, iuR, ipR = 0.5, 1.0, 0.125\n\nrho_u_p_LR = [irhoL, iuL, ipL, irhoR, iuR, ipR]\n\nx_ic = x_grid[id_ic, 0][:, None]\n\nt_ic = t_grid[id_ic, 0][:, None]\n\nx_ic_train = np.hstack((t_ic, x_ic))\n\nU_ic = IC(x_ic_train, irhoL, iuL, ipL, irhoR, iuR, ipR)\n\n#tensor\n\nx_ic_train = torch.tensor(x_ic_train, dtype=dtype, requires_grad=True, device=device)\n\nU_ic = torch.tensor(U_ic, dtype=dtype, device=device)\n\n#loss_xic(self, t_ic, x_ic):\n\nt_ic_train = to_numpy(x_ic_train)\n\nx_ic = t_ic_train[:,0:1]*1\n\n#tensor\n\nt_ic_train = torch.tensor(t_ic_train, dtype=dtype, device=device, requires_grad=True)\n\nx_ic = torch.tensor(x_ic, dtype=dtype, device=device)",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:44.631526Z",
     "iopub.execute_input": "2024-10-22T08:06:44.632252Z",
     "iopub.status.idle": "2024-10-22T08:06:44.643136Z",
     "shell.execute_reply.started": "2024-10-22T08:06:44.632211Z",
     "shell.execute_reply": "2024-10-22T08:06:44.642208Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:42.309631Z",
     "start_time": "2024-10-22T06:51:42.292094Z"
    }
   },
   "execution_count": 121,
   "outputs": []
  },
  {
   "id": "b5ad1a076b1101fa",
   "cell_type": "code",
   "source": "#残差\n#loss_pde(self, x_int):\nnum_f_train = 5000\nid_f = np.random.choice(num_x * num_t, num_f_train, replace=False)\nx_int = X[:, 0][id_f, None]\nt_int = T[:, 0][id_f, None]\nx_int_train = np.hstack((t_int, x_int))\n#tensor\nx_int_train = torch.tensor(x_int_train, dtype=dtype, requires_grad=True, device=device)",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:45.105374Z",
     "iopub.execute_input": "2024-10-22T08:06:45.105749Z",
     "iopub.status.idle": "2024-10-22T08:06:45.115341Z",
     "shell.execute_reply.started": "2024-10-22T08:06:45.105713Z",
     "shell.execute_reply": "2024-10-22T08:06:45.114405Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:42.874944Z",
     "start_time": "2024-10-22T06:51:42.863756Z"
    }
   },
   "execution_count": 122,
   "outputs": []
  },
  {
   "id": "b889c0c6fa2ecdf9",
   "cell_type": "code",
   "source": "#rh关系\n\n#loss_rh(self, t, x_start, x_end)\n\nt_rh_int_train = to_numpy(x_int_train)\n\nx_start = Xstart\n\nx_end = Xend\n\n#tensor\n\nt_rh_int_train = torch.tensor(t_rh_int_train, dtype=dtype, requires_grad=True, device=device)\n\nx_start = torch.tensor(x_start, dtype=dtype, device=device)\n\nx_end = torch.tensor(x_end, dtype=dtype, device=device)",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:45.338372Z",
     "iopub.execute_input": "2024-10-22T08:06:45.339078Z",
     "iopub.status.idle": "2024-10-22T08:06:45.344936Z",
     "shell.execute_reply.started": "2024-10-22T08:06:45.339042Z",
     "shell.execute_reply": "2024-10-22T08:06:45.343959Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:43.576092Z",
     "start_time": "2024-10-22T06:51:43.558603Z"
    }
   },
   "execution_count": 123,
   "outputs": []
  },
  {
   "id": "70089af1370df782",
   "cell_type": "code",
   "source": "#实例化神经网络\n\nlayers1 = Dnnlayers(2, 3, 3, 30)\n\nlayers2 = Dnnlayers(2, 3, 3, 30)\n\nlayers3 = Dnnlayers(2, 1, 1, 10)\n\n# Pinns(layers)\n\nmodel = conbine(layers1, layers2, layers3)\n\nloss_total_history = []\n\nloss_pde_history = []\n\nloss_ibc_history = []\n\nloss_rh_history = []\n\nlr = 0.001\n\noptimizer = torch.optim.Adam(list(model.model1.dnn.parameters())+list(model.model2.dnn.parameters())+list(model.model3.parameters()), lr=lr)",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:45.851854Z",
     "iopub.execute_input": "2024-10-22T08:06:45.852182Z",
     "iopub.status.idle": "2024-10-22T08:06:45.940322Z",
     "shell.execute_reply.started": "2024-10-22T08:06:45.852150Z",
     "shell.execute_reply": "2024-10-22T08:06:45.939550Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:44.217453Z",
     "start_time": "2024-10-22T06:51:44.207428Z"
    }
   },
   "execution_count": 124,
   "outputs": []
  },
  {
   "id": "f674917808bcff0d",
   "cell_type": "code",
   "source": "#training\n\ndef train(epoch):\n\n    model.model1.dnn.train()\n\n    model.model2.dnn.train()\n\n    model.model3.train()\n\n    def closure():\n\n        optimizer.zero_grad()     \n\n        \n\n        loss_pde = model.loss_pde(x_int_train)\n\n        loss_ic = model.loss_ic(x_ic_train, U_ic)\n\n        loss_xic = model.loss_xic(x_ic)\n\n        loss_rh = model.loss_rh(t_rh_int_train, x_start, x_end)\n\n        #loss_x = model.loss_x()\n\n                   \n\n        loss = loss_pde + 10*(loss_ic + loss_rh + loss_xic)#  + loss_x\n\n        if epoch%100==0:\n\n            print(f'epoch:{epoch}, loss:{loss:.8f},loss_pde:{loss_pde:.8f}, loss_IC:{loss_ic:.8f}, loss_xic:{loss_xic:.8f}, loss_rh:{loss_rh:.8f}')\n\n        if epoch%1000==0:\n\n            t_rh = t_rh_int_train[:, 0:1]\n\n            x_exact = t_rh*1\n\n            plt.figure()\n\n            plt.plot(to_numpy(x_predict), to_numpy(t_rh), 'r', label='x_predict')\n\n            plt.plot(to_numpy(x_exact), to_numpy(t_rh), 'b', label='x_exact')\n\n            plt.xlabel('x')\n\n            plt.ylabel('t')\n\n            plt.legend()\n\n            plt.show()\n        \n\n        loss_total_history.append(to_numpy(loss))\n\n        loss_pde_history.append(to_numpy(loss_pde))\n\n        loss_ibc_history.append(to_numpy(loss_ic+loss_xic))\n\n        loss_rh_history.append(to_numpy(loss_rh))\n\n        \n\n        loss.backward()\n\n        return loss\n\n    loss = optimizer.step(closure)",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:46.355610Z",
     "iopub.execute_input": "2024-10-22T08:06:46.356006Z",
     "iopub.status.idle": "2024-10-22T08:06:46.366895Z",
     "shell.execute_reply.started": "2024-10-22T08:06:46.355969Z",
     "shell.execute_reply": "2024-10-22T08:06:46.365966Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:44.840362Z",
     "start_time": "2024-10-22T06:51:44.823698Z"
    }
   },
   "execution_count": 125,
   "outputs": []
  },
  {
   "id": "59e6753e762b2257",
   "cell_type": "code",
   "source": "epochs1 = 15000\n\n#for training\n\ntic = time.time()\n\nfor epoch in range(1, epochs1+1):\n\n    #初值\n\n    x_predict_ic = (model.model3.forword(t_ic_train))\n\n    x_predict = model.model3.forword(t_rh_int_train)\n\n    #训练\n\n    train(epoch)    \n\n    #更新间断出入\n\n    t_ic_train = to_numpy(torch.hstack((t_ic_train[:,0:1], x_predict_ic[:,0:1])))\n\n    t_ic_train = torch.tensor(t_ic_train, dtype=dtype, device=device, requires_grad=True)\n\n    t_rh_int_train = to_numpy(torch.hstack((t_rh_int_train[:,0:1], x_predict[:,0:1])))\n\n    t_rh_int_train = torch.tensor(t_rh_int_train, dtype=dtype, device=device, requires_grad=True)\n\ntoc = time.time()",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:06:46.876401Z",
     "iopub.execute_input": "2024-10-22T08:06:46.876751Z",
     "iopub.status.idle": "2024-10-22T08:14:34.595572Z",
     "shell.execute_reply.started": "2024-10-22T08:06:46.876715Z",
     "shell.execute_reply": "2024-10-22T08:14:34.594652Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:45.592121Z",
     "start_time": "2024-10-22T06:51:45.384310Z"
    }
   },
   "execution_count": 126,
   "outputs": []
  },
  {
   "id": "cdd6739e70b2ed04",
   "cell_type": "code",
   "source": "optimizer = torch.optim.LBFGS(list(model.model1.dnn.parameters())+list(model.model2.dnn.parameters())+list(model.model3.parameters()),lr=1.0,max_iter=1)\n\nepochs2 = 500\n\ntic = time.time()\n\nfor epoch in range(1, epochs2+1):\n\n    #初值\n\n    x_predict_ic = (model.model3.forword(t_ic_train))\n\n    x_predict = model.model3.forword(t_rh_int_train)\n\n    #训练\n\n    train(epoch)\n\n    #更新间断出入\n\n    t_ic_train = to_numpy(torch.hstack((t_ic_train[:,0:1], x_predict_ic[:,0:1])))\n\n    t_ic_train = torch.tensor(t_ic_train, dtype=dtype, device=device, requires_grad=True)\n\n    t_rh_int_train = to_numpy(torch.hstack((t_rh_int_train[:,0:1], x_predict[:,0:1])))\n\n    t_rh_int_train = torch.tensor(t_rh_int_train, dtype=dtype, device=device, requires_grad=True)\n\n    \n\ntoc = time.time()",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:18:33.267417Z",
     "iopub.execute_input": "2024-10-22T08:18:33.268029Z",
     "iopub.status.idle": "2024-10-22T08:18:49.434546Z",
     "shell.execute_reply.started": "2024-10-22T08:18:33.267987Z",
     "shell.execute_reply": "2024-10-22T08:18:49.433481Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:46.991583Z",
     "start_time": "2024-10-22T06:51:46.794825Z"
    }
   },
   "execution_count": 132,
   "outputs": []
  },
  {
   "id": "6f74f9ec3398306c",
   "cell_type": "code",
   "source": "#test\nx_predict = (model.model3.forword(t_rh_int_train))\ndx = gradients(x_predict, t_rh_int_train)[0]\ndx_t = dx[:,0:1]\nprint(dx_t)\ndx_x = dx[:,1:2]\nprint(dx_x)\nprint(x_predict-t_rh_int_train[:,1:2])\nx_predict = to_numpy(model.model3.forword(t_rh_int_train))\n\nt_rh = to_numpy(t_rh_int_train[:,0:1])\n\nx_exact = t_rh*1.0\n\nplt.figure()\n\nplt.plot(x_predict, t_rh, 'r', label='x_predict')\n\nplt.plot(x_exact, t_rh, 'b', label='x_exact')\n\nplt.xlabel('x')\n\nplt.ylabel('t')\n\nplt.legend()\n\nplt.show()\n\nerror = abs(x_predict - x_exact)\n\nplt.figure()\n\n#plt.plot(t_rh, error, label='x_error')\n\nplt.scatter(t_rh, error, label='x_error', s=0.5)\n\nplt.xlabel('t')\n\nplt.ylabel('error')\n\nplt.legend()\n\nplt.show()",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:18:54.956430Z",
     "iopub.execute_input": "2024-10-22T08:18:54.956806Z",
     "iopub.status.idle": "2024-10-22T08:18:55.486654Z",
     "shell.execute_reply.started": "2024-10-22T08:18:54.956770Z",
     "shell.execute_reply": "2024-10-22T08:18:55.485734Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:48.898191Z",
     "start_time": "2024-10-22T06:51:48.580874Z"
    }
   },
   "execution_count": 133,
   "outputs": []
  },
  {
   "id": "5dae7274de0cce36",
   "cell_type": "code",
   "source": "x_int_test = to_numpy(x_int_train)*1.0\n\n#exact\n\nrho_exact = np.zeros_like(x_int_test[:, 0:1])\n\nu_exact = np.zeros_like(x_int_test[:, 0:1])\n\np_exact = np.zeros_like(x_int_test[:, 0:1])\n\nu_exact[:, 0:1] = 1.0\n\np_exact[:, 0:1] = 0.125\n\n\n\n#x_int_test = torch.tensor(x_int_test, dtype=dtype, requires_grad=True, device=device)\n\nx = x_int_test[:, 1]\n\n\n\nindex1 = np.where(x <= x_predict[:,0])\n\nindex2 = np.where(x > x_predict[:,0])\n\nrho_exact[index1] = 1.0\n\nrho_exact[index2] = 0.5\n\n\n\nx1_int = x_int_test[index1]\n\nx2_int = x_int_test[index2]\n\nx1_int = torch.tensor(x1_int, dtype=dtype, device=device, requires_grad=True)\n\nx2_int = torch.tensor(x2_int, dtype=dtype, device=device, requires_grad=True)\n\n\n\nrho_predict1, u_predict1, p_predict1 = (model.model1.forword(x1_int))\n\nrho_predict2, u_predict2, p_predict2 = (model.model2.forword(x2_int))\n\n\n\nrho_predict = to_numpy(torch.vstack([rho_predict1, rho_predict2]))\n\nrho_predict[index1] = to_numpy(rho_predict1)\n\nrho_predict[index2] = to_numpy(rho_predict2)\n\nu_predict = to_numpy(torch.vstack([u_predict1, u_predict2]))\n\nu_predict[index1] = to_numpy(u_predict1)\n\nu_predict[index2] = to_numpy(u_predict2)\n\np_predict = to_numpy(torch.vstack([p_predict1, p_predict2]))\n\np_predict[index1] = to_numpy(p_predict1)\n\np_predict[index2] = to_numpy(p_predict2)\n\n\n\n'''t_grid_test = to_numpy(x_int_test[:, 0:1]).reshape(num_x, num_t)\n\nx_grid_test = to_numpy(x_int_test[:, 1:2]).reshape(num_x, num_t)\n\nrho_predict = to_numpy(rho_predict).reshape(num_x, num_t)\n\nu_predict = to_numpy(u_predict).reshape(num_x, num_t)\n\np_predict = to_numpy(p_predict).reshape(num_x, num_t)\n\nrho_exact = rho_exact.reshape(num_x, num_t)\n\nu_exact = u_exact.reshape(num_x, num_t)\n\np_exact = p_exact.reshape(num_x, num_t)'''\n\n\n\nerror_rho = abs(rho_predict - rho_exact)\n\nerror_u = abs(u_predict - u_exact)\n\nerror_p = abs(p_predict - p_exact)\n\n# 绘制带有颜色映射的散点图\n\nplt.figure()\n\nplt.scatter(x_int_test[:, 1:2], x_int_test[:, 0:1], c=rho_predict, cmap='viridis', alpha=0.8)\n\nplt.colorbar()\n\nplt.xlabel('x')\n\nplt.ylabel('t')\n\nplt.title('rho_predict')\n\nplt.show()\n\nplt.figure()\n\nplt.scatter(x_int_test[:, 1:2], x_int_test[:, 0:1], c=u_predict, cmap='viridis', alpha=0.8)\n\nplt.colorbar()\n\nplt.xlabel('x')\n\nplt.ylabel('t')\n\nplt.title('u_predict')\n\nplt.show()\n\nplt.figure()\n\nplt.scatter(x_int_test[:, 1:2], x_int_test[:, 0:1], c=p_predict, cmap='viridis', alpha=0.8)\n\nplt.colorbar()\n\nplt.xlabel('x')\n\nplt.ylabel('t')\n\nplt.title('p_predict')\n\nplt.show()\n\nplt.figure()\n\nplt.scatter(x_int_test[:, 1:2], x_int_test[:, 0:1], c=error_rho, cmap='viridis', alpha=0.8)\n\nplt.colorbar()\n\nplt.xlabel('x')\n\nplt.ylabel('t')\n\nplt.title('error_rho')\n\nplt.show()\n\nplt.figure()\n\nplt.scatter(x_int_test[:, 1:2], x_int_test[:, 0:1], c=error_u, cmap='viridis', alpha=0.8)\n\nplt.colorbar()\n\nplt.xlabel('x')\n\nplt.ylabel('t')\n\nplt.title('error_u')\n\nplt.show()\n\nplt.figure()\n\nplt.scatter(x_int_test[:, 1:2], x_int_test[:, 0:1], c=error_p, cmap='viridis', alpha=0.8)\n\nplt.colorbar()\n\nplt.xlabel('x')\n\nplt.ylabel('t')\n\nplt.title('error_p')\n\nplt.show()\n\n\n\n\n\n\n\nprint('L2 for error_rho is :', np.sqrt(((error_rho)**2).mean()))\n\nprint('L2 for error_u is :', np.sqrt(((error_u)**2).mean()))\n\nprint('L2 for error_p is :', np.sqrt(((error_p)**2).mean()))\n\n'''\nL2 for error_rho is : 0.0005062342402445959\nL2 for error_u is : 0.0006041250081679045\nL2 for error_p is : 0.00035317003769470416\n'''",
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-10-22T08:18:55.488684Z",
     "iopub.execute_input": "2024-10-22T08:18:55.489118Z",
     "iopub.status.idle": "2024-10-22T08:18:58.660480Z",
     "shell.execute_reply.started": "2024-10-22T08:18:55.489072Z",
     "shell.execute_reply": "2024-10-22T08:18:58.659556Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-10-22T06:51:51.330357Z",
     "start_time": "2024-10-22T06:51:49.981717Z"
    }
   },
   "execution_count": 134,
   "outputs": []
  },
  {
   "id": "f56e33350875c610",
   "cell_type": "code",
   "source": "",
   "metadata": {},
   "execution_count": null,
   "outputs": []
  }
 ]
}
